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Abstract 

The generalized linear model is widely used in all areas of applied statistics 
and while correct asymptotic inference can be achieved under misspecihca- 
tion of the distributional assumptions, a correctly specihed mean structure 
is crucial to obtain interpretable results. Usually the linearity and func¬ 
tional form of predictors are checked by inspecting various scatterplots of 
the residuals, however, the subjective task of judging these can be challeng¬ 
ing. In this paper we present an implementation of model diagnostics for 
the generalized linear model as well as structural equation models, based on 
aggregates of the residuals where the asymptotic behavior under the null is 
imitated by simulations. A procedure for checking the proportional hazard 
assumption in the Cox regression is also implemented. 
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1. Introduction 


The generalized linear model is one of the most widely used classes of 
statistical models, however, the standard methods of inference relies on dis¬ 
tributional and linearity assumptions. The importance of this is sometimes 
underestimated, to some extent because few tools are available for checking 
all the aspects of the model. While the distributional assumptions can be 
relaxed, i.e., by using a sandwich estimator as implemented in the sandwich 


package (Zeileis, 2006), careful attention should be paid to the validity of 
the specihed mean structure. A typical model check involves assessment 
of various residual plots. As the true variance of individual residuals are 
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unknown it can be difficult to decide whether a residual plot indicates a 
reasonable specihcation of the mean or not. In a paper by |Su and Wei 


(1991) it was proposed instead to look at certain aggregates of the resid¬ 
uals, such as the cumulative sum over predicted values or covariates. The 
key result here is, that the asymptotic distribution of such aggregates can 
be determined under the hypothesis that the model is correctly specihed. 


The R environment is one of the most widely used statistics platforms 
but lacks objective diagnostics tools for many regression models, and in 
particular methods based on aggregates of residuals, thus motivating the 
creation of the gof-package described in the following sections. 


2. Implementation 

The gof package implements diagnostics of the linearity assumptions 
for the generalized linear model and linear structural equation models. Fur¬ 
ther similar methods are available for checking the proportional hazards 
assumption of the Cox regression model for right censored data. The fol¬ 
lowing section describes the theoretical details behind the implementation. 


2.1. Generalized linear model 


The case of generalized linear models was hrst examined by Su and Wei 


(1991). Let Y be the response variable with a distribution from a (natural) 
exponential family: 


f{Y = Vi I 0) = exp I + c{yi, 0) 


( 1 ) 


parameterized by 9 (and the dispersion parameter 0) and the known func¬ 
tions a, b and c. Direct calculations reveals that the 


EDi = b'{ei), Var(y,) = a(0)6"(0,). (2) 

The mean El^ = is related to some covariates, Xi, through a link- 

function (McCullagh and Nelder, 1983), g, 

9 = 0 ^ Xi, (3) 


i.e., 9i = 6i{f3). Typically, the canonical link is chosen such that o p = id, 
with the most common regression models being the general linear model, 
logistic regression and Poisson regression 
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family canonical link a(0) b'{6) 

Normal identity 0 9 

Binomial logit 1 1/(1 + exp(—0)) 

Poisson log 1 exp (6*) 

Given n observations (?/,, Xu ,..., the maximnm likelihood es¬ 

timate f3 eW is obtained by solving the set of score eqnations: 


U{(3) = '^h{f3^Xi)xi{yi-g \f3^Xi)] , (4) 

i=l 


with h = d{{g o y) We dehne the (raw) residnals Ci = yi — g ^{fB'^Xi), 
i = 1,... ,n. Onr interest is the cnmnlative snm of the residnals over the 


jth covariate (Sn and Wei, 1991 Lin et ah, 2002): 


Wj{x) =n ^ (5) 

i=l 


In contrast to the distribntion of individnal residnals, we can determine 
the distribntion (nnder the nnll) of this aggregate. For known parameters 
the a symptotics can be derived as a Brownian bridge (|Shorack and Wellner 


1986), however, we need to take nncertainty in estimation of (3 into acconnt. 


Under certain regnlarity conditions, a Taylor expansion aronnd the trne 
parameter valne, /3o, gives ns 


W^{x) = W,{x I 3) = W,{x I /3o) + I /3) 


(/3 — /3o) + Op(l). 


/3=/3o 


( 6 ) 


Let X(/3) = E(—Vt/(/3)) denote the Fisher information. Now (/3 — /3o) 
is asymptotically normally distribnted and asymptotically eqnivalent with 

i0)-^u0y. 

0-(3o)=AP)-^U0)+o,{1). (7) 


It then follows that the process 


Wj{x) =n ^ + vM I ^{fi)Xih{p'^Xi) 

i=l 


i 


( 8 ) 
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with i.i.d. Gi,..., Gn ~ A/'(0,1), i = 1,..., n, and 


Vjix I (3) 


n 

2 = 1 


df3 


( 9 ) 


(see Table converges weakly to the same limiting distribution as the 
observed process ([^ (Lin et al., 2002). 


9{x) 

9 \z) 

d{9-^){z) 

X 

z 

1 

logit (x) 

l/(l + exp(-z)) 

exp(—^)/[l -|- exp(—^)]^ 

log(x) 

exp(z) 

exp(z) 


Table 1: Some link functions and their inverse. 

To test the functional form of the jth covariate we look at a Kolmogorov- 
Smirnov (KS) type supremum statistic: 

‘3:^£-Wj^snp\Wj{x)\. ( 10 ) 

X 

Alternatively tests can be based on the Cramer-von-Mises (CvM) functional: 

‘^ 2 ^- ^ [\Wj{x)f dx. (11) 


A large number of realizations of Wj is generated. The supremum statistic 
is calculated for each realization and the p-value is estimated from the em¬ 
pirical distribution of these statistics. The residuals can also be cumulated 


after the predicted values (Lin et ah, 2002) 




(l3^Xi)<t} 


^2? 


( 12 ) 


2=1 


which leads to a test of misspecihed link function. 


2.2. Structural equation models 

The linear structural equation models covers a broad range of models 
including the general linear model, path analysis and various latent variable 
models. Diagnostics based on cumulative residuals was examined in this 


case by Sanchez et al. (2009) building on the work of Pan and Lin (2005) on 
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Generalized Linear Mixed Models (GLMM) sharing many of the aspects of 
structural equation models. The basic idea and proof of weak convergence 
is very similar to the case of GLM. 

A structural equation model is typically divided into two separate parts. 
For the Ah individual we have a measurement part describing the multivari¬ 
ate outcome Yp 


Yi — jy + Ar]i + KXi + Si, (13) 

where rfi are the latent variables and Xi are covariates, and a structural 
part describing the latent variables: 


rji — ex + Brji + TXi -|- (14) 

where v G MP, A G K G and e* ~ Ap(0, Eg). And ex G 

B G r G and Ci ~ Ar(0,^'). Hence, the model is parameter¬ 

ized by some 6 dehning (iz, ex, A, K, B, F, Eg, ^') with some restrictions to 
guarantee identiheation. The conditional moments of T) given X^, are 

^i = E 0 {Yi\ Xi) = IX + A{1 - B)~^ex 

+ [A(l - B)-^T + K] Yi, 

E = Vare(r, I X,) = A(1 - B)-^^{1 - B)-^^A^, (16) 

and inference on 6 is usually obtained by MLE [Bollen 

The residuals can be predicted as the conditional mean given the en¬ 
dogenous variables and covariates. Hence, 

Cik = E{eik I Yi, Xi) = 7rfEgE~^(T) - Hi) (17) 

5, = E(C, I Y„ X,) = - B)-^^A^'E-\Yi - pa) (18) 

where W —)■ M is the projection onto coordinate s. Different local 
aspects of the structural equation model can now be assessed by examining 
the cumulative residual processes of either 'cik and Qg. 

Misspecified covariate effect on the gth latent variable is checked by sum¬ 
ming Qg with respect to the jth covariate, {Xij): 

n 

^ ( 19 ) 


(1989). 
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and as for the GLM we can imitate the behavior of this process under the 
null of no misspecihcation by simulation. Misspecified link between gth latent 
variable and its predictors is checked by summing (^ig with respect to \ 
Xi) = 77^(1 — + 'yXi). To examine departures from the specified 

association between an endogenous variable and one of its predictors, we 
can look at the cumulative process dehned by summing Zik with respect to 
Xj or E,{pig I Xi). This can also be used to diagnose for so-called item 
bias (conditional dependence between a covariate and endogenous variable 
given latent variables). Finally, misspecified link between an endogenous 
variable and its linear predictors is checked by summing with respect to 

ny^k I ^*). 


2.3. Cox’s proportional hazard model 

The idea of looking at aggregates of residuals can also be applied as 
a tool for diagnosing the proportional hazards assumptions used in many 
survival analyses. We will assume that we have triplet observations {Ni{t), 
Yi{t), Xi{t)), i = 1,..., n of a counting process, at-risk process and covariate 
process in the compact time-interval [0,r]. Using the notation of stochastic 
integrals we let the Martingale decomposition of the counting process be 
given by 

dNi{t) = Xi{t) dt + dMiit). (20) 


Cox’s proportional hazard model assumes intensity takes the form 

Xiit) = U(t)Ao(t)exp(Xf(f)/?), (21) 

where X is p-dimensional covariates. We denote the cumulative baseline 
hazard 


Ao(i) ^ / \o{s)ds. 


( 22 ) 


As the model contains a non-parametric term, Aq, inference will be based 
on the partial likelihood (Cox, 1972[ ) 

p3) 


i=l t 


where 


So{t,^) = ^U(t)exp(Xf(f)/3). 


(24) 
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with the hrst and second partial derivatives 


= ^yi(t)exp(Xf(t)/3)Xi(t), 

i 

S2{t,l3) = X^r.(«)exp(Jff(i)«V(i)®t 

i 

and let E(t,l3) = Si/So(t, /3). The score equation then becomes 
UW = ^ [ |A',(i) - E(t,P)] dN,(t). 


(25) 

( 26 ) 


( 27 ) 


i.i •'» 


The Nelson-Aalen estimator of the cumulative intensity is 


Ao(t) — 


/o So{s,/3) 


dNXs), 


( 28 ) 


where iV. = Dehne 

71 rt 


i=i do 


/^(s,li)-E{s,Pf^dN,(s)= [ V{s,li)dNXs), ( 29 ) 

*^0 Jo 


and hence minus the derivative of the score is /(r, P) (i.e., the information). 
The estimated martingales residual process is given by 


Mi{t) = Ni{t) - Ai{t) = Ni{t) - / Fi(s) exp Xl{s)p] dAo(s) 


= Ni{t) - / exp Xf {s)l3 


(30) 


So{s,P) 


dNXs), 


(with the martingale residuals dehned by evaluation in r), and the estimated 
score process 


71 rt 


uip,t) = y 2 / xxs)dMps) 

i=l do 

n 

= E/ {x,M-E(s,^)] dN,(s), 


( 31 ) 


i=l -^0 


where Xi(s) — E{s,P) are the Sehoenfeld residuals. 
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To assess the proportional hazards assumption we will calculate the 
Kolmogorov-Smirnov and Cramer-von-Mises test statistics of the different 
coordinates of the observed score process. As in the previous section we can 
simulate realizations under the null (proportional hazards). The key result 
is that is asymptotically equivalent to 

OO 

(32) 

i=l 

with 


Muit)= / {V(s)-e(s,/?o)}dMi(a), 


(33) 


P 

where e{t,(3o) = limn^ao E{t, (3 q) (see (Martinussen and Scheike 


2006 


Lin 


et al.f 1993)), which follows from a Taylor expansion around the true pa¬ 


rameter (3q. With the estimates plugged in we get 

Mu{t)= [ {Xi{s)-E{sJ))dN,{s) 

Jo 

exp(Xf(s)/3) 


(W(s)-E(s,/3))- 


SoisJ) 


dNXs). 


(34) 


Given observed times (Ti,... ,T„) and death-indicators Aj, (Xj = Xi(Tj)) 
we can implement this by 


l(T,<<,A=i) |.V(s) - E{T„ J)| - exp(XTfl)Mt) 

+ exp{X[D) f Si/S'i{s,P)dN,{s). 

Jo 

(36) 


Finally n is asymptotically equivalent to n Mii(t)Gi where 

the Gj’s are i.i.d. 7V(0,1). 


2.J. Software 

The described methods are implemented in the R-package gof available 


from the Comprehensive R Archive Network (R Core Team, 2012). 


The package has been designed to work directly on Im, glm and coxph 
objects (Therneau and original R port by Thomas Lumley, 2013). Addition¬ 
ally, various aspects of latent variable models, htted via the lava-package 
(Holst and Budtz-Joergensen| 2012), can be diagnosed. 


















The simulation routine is computational intensive and to obtain bet¬ 
ter computing efficiency, the resampling routines was written in C++. The 


implementation uses the Scythe Statistical Library (Pemstein et ah 


2011) which among other things offers operator overloaded matrix opera¬ 
tions making the (linear) algebraic computations in the program close to 
self-documenting. 


3. Examples 

In the following section the gof package will be demonstrated in gen¬ 
eralized linear models, a structural equation model and a Cox regression 
model. 

3.1. Generalized linear models 

First we define a simple function that allows us to simulate data from 
Binomial and Poisson regression models with link function and covariates 
X, Z ~ A^(0,1) 


g{¥.[Y\X,Z]) = f{X,Z). (36) 


R> siml <- function(n,f=sum,family=binomial("logit")) { 

X <- rnorm(n) 
z <- rnorm(n) 

if (is.character(family)) family <- do.call(family,list()) 
eta <- family$linkinv(apply(cbind(x,z),l>£)) 
y <- switch(family$family, 

binomial= (eta>runif(n))*1, 
poisson= rpois (n,eta), 
eta) 

return(data.frame(y,x,z)) 


} 


We first simulate binomially distributed observations and use a comple¬ 
mentary log-log link: 


log(-log[l-E(F |X,Z)]) =X + Z (37) 


R> d <- siml(n=1000,family=binomial("cloglog")) 
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Next we fit both the correct model and the model with canonical link 

R> 11 <- glm(y~x+z,d,family=binomial("cloglog")) 

R> 12 <- glm(y~x+z,d,family=binomial("logit")) 

Using the cumres method, we calculate the cumulative residual process 
ordered by the predicted values and simulate 1,000 processes from the null 

R> library("gof") 

R> (gl <- cumres(11,R=1000,variable="predicted")) 

Kolmogorov-Smirnov-test: p-value=0.466 
Cramer von Mises-test: p-value=0.333 

Based on 1000 realizations. Cumulated residuals ordered by predicted-variable 


Kolmogorov-Smirnov-test: p-value=0.466 
Cramer von Mises-test: p-value=0.333 

Based on 1000 realizations. Cumulated residuals ordered by predicted-variable 


R> (g2 <- cumres(12,R=1000,variable="predicted")) 

Kolmogorov-Smirnov-test: p-value=0.014 
Cramer von Mises-test: p-value=0 

Based on 1000 realizations. Cumulated residuals ordered by predicted-variable 


Kolmogorov-Smirnov-test: p-value=0.014 
Cramer von Mises-test: p-value=0 

Based on 1000 realizations. Cumulated residuals ordered by predicted-variable 


There are clear indications, by both the supremum and CvM test, of mis- 
specihcation of the link function in model 12. To plot the observed pro¬ 
cess and realizations from under the null (the number of realization can be 
changed in the cumres call with the argument plots), we can use the plot 
method 
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R> par(mfrow=c(l,2)) 

R> plot(gl,title="Model 'll"'); plot(g2,title="Model '12'") 


Model 'll' 


Model '12' 


KS-test: p.0.466 
CvM-lesI: p.0.333 




1-i-1-1-1-r 

0.0 0.2 0.4 0.6 0.8 1.0 

predicted 



Figure 1: Cumulative residual processes of model 11 and 12 with residuals ordered by 
the predicted response. The gray curves are 50 realizations from the null model. The 
transparent blue area defines a 95% prediction band for all the simulated processes. 

It is evident from the plot (Figure ??), that the observed process of 
model 2 is extreme. 

Next we simulate data from a Poisson regression model 

\og{E{Y \ X,Z)) =0.5-+ Z (38) 


R> d2 <- siml(200,f=f unction (x) 0.5*x[l]~2+x[2] ,family=poisson()) 

and we £t a Poisson regression model but with misspecihed functional form 
of the covariate X 

R> 1 <- glm(y~x+z,family=poisson(),data=d2) 

Next we check the link function and functional form of both covariates 

R> (g <- cumres(l,R=2000)) 

Kolmogorov-Smirnov-test: p-value=0.453 
Cramer von Mises-test: p-value=0.547 

Based on 2000 realizations. Cumulated residuals ordered by predicted-variable. 

Kolmogorov-Smirnov-test: p-value=0.001 
Cramer von Mises-test: p-value=0 

11 







Based on 2000 realizations. Cumulated residuals ordered by x-variable. 

Kolmogorov-Smirnov-test: p-value=0.8185 
Cramer von Mises-test: p-value=0.858 

Based on 2000 realizations. Cumulated residuals ordered by z-variable. 


Kolmogorov-Smirnov-test: p-value=0.453 
Cramer von Mises-test: p-value=0.547 

Based on 2000 realizations. Cumulated residuals ordered by predicted-variable. 

Kolmogorov-Smirnov-test: p-value=0.001 
Cramer von Mises-test: p-value=0 

Based on 2000 realizations. Cumulated residuals ordered by x-variable. 

Kolmogorov-Smirnov-test: p-value=0.8185 
Cramer von Mises-test: p-value=0.858 

Based on 2000 realizations. Cumulated residuals ordered by z-variable. 


and we plot all processes (Figure while changing the color (and alpha 
blending) of the realizations and prediction-band (setting col or col.ci to 
NULL will disable either the realizations or the prediction-band) 
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R> par(mfrow=c(2,2)) 

R> plot (g,col="gray",col.ci = "black",col.alpha=0.4,legend=NULL) 



X 

s' 



z 


Figure 2: Cumulative residual processes for the Poisson regression model 1. 


Again, the misspecification (of the functional form of X) is evident from 
the plots. 

3.2. Structural equation models 

The cumres method is also available for structural equation models htted 


via the lava package (Holst and Budtz-Joergensen, 2012). As an example 


we will examine a simple model, with three outcomes described by the 
equation 


Pj T T 6jj, j 1,..., 3, 


(39) 
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with i = 1,..., n individuals and latent variable rji. We also add a structural 
equation describing the latent variable 


V^ = Pl■X + ^2Z + C, (40) 

with covariates X and Zl. The residual terms en,... ,€i 3 X cire normally 
distributed and independent. In lava we can specify the model as 

R> library(lava) 

R> m <- Ivm(list(c(yl,y2,y3)~eta,eta~x+z)) 

R> latent (m) <- ~eta 

We simulate 200 observations from a structural equation model like the one 
dehned above, with intercepts set to zero and all other parameters equal to 
one, but with 


Yi 2 — 1 ]“^ €i 2 and — X 0.5 ■ X"^ Z -j- (^. (41) 


R> mO <- m 

R> functional(m0,y2~eta) <- function(x) x~2 
R> functional(mO,eta~z) <- function(x) x+0.5*x~2 
R> d <- sim(mO,200) 

Next we hud the MLE of the hrst model 

R> (e <- estimate(ni,d)) 



Estimate 

Std. Error 

Z-value 

P-value 

Measurements: 

y2<-eta 

2.21830 

0.24535 

9.04147 

<le-12 

y3<-eta 

0.99314 

0.05380 

18.45847 

<le-12 

Regressions: 

eta<-x 

1.00280 

0.09820 

10.21176 

<le-12 

eta<-z 

1.08055 

0.09790 

11.03729 

<le-12 

Intercepts: 

y2 

2.84517 

0.48757 

5.83545 

5.365e-09 

y3 

-0.05982 

0.09412 

-0.63558 

0.5251 

eta 

0.50307 

0.10623 

4.73553 

2.185e-06 

Residual Variances: 

yi 

0.67836 

0.15287 

4.43739 


y2 

40.57804 

4.22914 

9.59488 


y3 

0.92828 

0.16468 

5.63682 


eta 

1.57358 

0.21605 

14 

7.28351 




and as an example we cumulate the predicted residual terms of and Y 2 
against E,{r]i \ Xj), and the residual term of T]i against the two covariates. 

R> e.gof <- cumres(e,list(y3~eta,y2~eta,eta~x,eta~z),R=1000) 

From the cumulative residual plots (see Figure we clearly see the mis- 
specihcation in the measurement model of the second outcome (with the 
observed process also indicating a quadratic form), and also the wrongly 
specihed functional form of X. 

For complete flexibility the cumres method can be used with the syntax 
cumres(model,y,x, . . .), where y is a function of the model parameters 
returning the residuals of interest, and x can be any vector to order the 
residuals by. Typically y will be dehned via the predict method of a 
Ivmf it object (a lava model object). 
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y3 


y2 


tn 

d 


KS-lest: p«0.7 
CvM-lest p-0.723 

8 - 

0 _ 


KS-test: p.0.003 
CvM-tesI: p-0.002 

X 

? ■ 


X 

0 - 

0 

CM - 
1 




-2 0 2 4 -2 0 2 4 


eta 


eta 


eta 


eta 




X 


Figure 3: par(mfrow=c( 2 , 2 )); lapplyCe.gof ,plot). Selected cumulative residual 
processes for the structural equation model fit e. The top row shows the cumulative 
residuals of £ 3 ^ and £ 2 ^ (see (39)) ordered by £( 77 ^ | Xi,Zi). The bottom row shows the 
cumulative processes of the predicted residual term, Q, of the latent variable ordered by 
each of the two covariates. 


3.3. Cox regression - Mayo clinic PBC data 

As an example of checking the proportional hazards assnmption in a Cox 
model, we will analyze the Mayo Clinic PBC data. Dickson et al. (1989) 


suggested a Cox model for analyzing the survival of the liver disease patient 
with 5 covariates: age, edema status, logarithmic serum bilirubin, logarith¬ 
mic standardized blood clotting time, and logarithmic serum albumin: 


R> library("survival") 

R> dataC'pbc") 

R> pbc.cox <- coxph(Surv(time,status==2)~age+edema+log(bili)+ 
log(protiine)+log(albumin), data=pbc) 
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To check the proportional hazards assumption, we examine the score 
process vs. follow-up time: 

R> pbc.gof <- cumres(pbc.cox,R=2000) 

and plot the observed process with realizations from the null 

R> par(mfrow=c(2,3)) 

R> plot(pbc.gof,legend=FALSE) 



age 


edema 


log(bili) 

8- 

8 - 

M 

■S o - 
5? 

8 _ 

1 

O 

a - 


—1-1 1-1 1— 

01 S 0 9- 01- 


~20 -10 0 10 20 

_1_1_1_1_1_ 



T—1—1—1—H 

0 1X0 3X0 

Tm« 

log( protime) 


T—1—1—1—r 

0 1X0 3X0 

Tm# 

log(albumin) 


T—1—1—1—H 

0 1X0 3X0 

Tm« 

M 

1 

B O - 

a 

9 

? CJ 

• 

1 


-1 1-1-1 1- 

Z lOl-Z- 

Vv^ 




T-1-1-1—H 

0 1X0 3X0 

Tni« 


1 1 1 1 

0 1X0 3X0 

Tnt« 




Figure 4: Cumulative score processes for the Cox regression analysis of the PBC data, 
pbc.cox. 

There are clear indication of violation of the proportional hazards as¬ 
sumption for blood clotting time (protime), and indication of problems with 
the edema variable. To remedy the non-proportionality, time-varying covari¬ 
ate effects could be introduced to the model, e.g., 

R> library("timereg") 

R> pbc.caalen <- cox.aalen(Surv(time,status==2) ~ prop(age) + prop(edema) 
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prop(bili) + protime, data=pbc, n.sim=500) 


4. Conclusion 

The package gof adds a valuable tool to the model diagnostics toolbox 
and gives an objective method for evaluating the linearity assumptions in the 
generalized linear model and linear structural equation models. Extensions 
to other models such as the linear mixed model can be implemented using 
the C++ interface as used by the cumres method for glm and Ivm objects. 
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